## estimate MQ-scores

library(MCMCpack)

## load data
load("mqData2013.Rda")
votes <- mqdata[,1:45]

## remove unanimous cases
indic <- rep(TRUE, nrow(votes))
for (i in 1:nrow(votes)){
  meanval <- apply(votes[i,], 1, mean, na.rm=TRUE)
  if (meanval == 0 | meanval == 1){
    indic[i] <- FALSE
  }
}
votes <- votes[indic,]
mqdata <- mqdata[indic,]
IT <- mqdata$time


## set tau
tau2.start <- rep(0.1, 45)

## set priors
e0 <- rep(0, 45)
E0 <- rep(1, 45)
e0[2] <- -2 ## Black
E0[2] <- 10
e0[4] <- 1 ## Stewart
E0[4] <- 10
e0[21] <- 3 ## Rehnquist
E0[21] <- 10

## set starting values
theta.start <- rep(0, 45)
theta.start[1] <- 1
theta.start[2] <- -2
theta.start[3] <- -2
theta.start[4] <- 1
theta.start[5] <- -1
theta.start[6] <- -1
theta.start[11] <- 1
theta.start[18] <- 1
theta.start[21] <- 3
theta.start[24] <- 3
theta.start[27] <- 3
theta.start[35] <- -1
theta.start[37] <- 1
theta.start[38] <- 1
theta.start[42] <- 1
theta.start[43] <- 1
theta.start[44] <- -1



  out <- MCMCdynamicIRT1d(t(votes),
                          item.time.map=IT,
                          theta.start=theta.start,
                          mcmc=1000000, burnin=25000, thin=100,
                          verbose=500, seed=24680, tau2.start=tau2.start,
                          e0=e0, E0=E0,
                          a0=0, A0=1,
                          b0=0, B0=1, c0=-10, d0=-2,
                          store.item=FALSE,
                          theta.constraints=list(Thomas="+", Alito="+",
                            Brennan="-", Douglas="-", Whittaker="+"))

save.image("MQout.Rdata")

